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Abstract — We present a pursuit-like algorithm tliat we 
call the "superset method" for recovery of sparse vectors 
from consecutive Fourier measurements in the super- 
resolution regime. The algorithm has a subspace iden- 
tification step that hinges on the translation invariance 
of the Fourier transform, followed by a removal step 
to estimate the solution's support. The superset method 
is always successful in the noiseless regime (unlike £i 
minimization) and generalizes to higher dimensions (unlike 
the matrix pencil method). Relative robustness to noise is 
demonstrated numerically. 
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I. Introduction 

We consider the problem of recovering a sparse vector 
xq G R", or an approximation thereof, from m < n 
contiguous Fourier measurements 

y = AxQ+e, (1) 

where A is the partial, short and wide Fourier matrix 
Ajk = 6^'^*^'=/", < j < m, -n/2 < fc < n/2, n even, 
and, say, e ~ iV(0, a'^Im)- 

When recovery is successful in this scenario of 
contiguous measurements, we may speak of super- 
resolution: the spacing between neighboring nonzero 
components in xq can be much smaller than the Rayleigh 
limit n/in suggested by Shannon-Nyquist theory. But in 
contrast to the compressed sensing scenario, where the 
m values of j are drawn at random from {0,...,ri— 1}, 
super-resolution can be arbitrarily ill-posed. Open ques- 
tions concern not only recovery bounds, but the very 
algorithms needed to define good estimators. 

Various techniques have been proposed in the liter- 
ature to tackle super-resolution, such as MUSIC ifTTI . 
Prony's method / finite rate of innovation ||8l H] lfT2l . 



the matrix pencil method Jg), £i minimization Q 15] ill 
Il2l, and greedy pursuits JU. 

Prony and matrix pencil methods are based on eigen- 
value computations: they work well with exact mea- 
surements, but their performance is poorly understood 
in the presence of noise, and they are not obviously 
set up in higher dimensions. As for £i minimization, 
there is good evidence that fc-sparse nonnegative signals 
can be recovered from only 2fc + 1 noiseless Fourier 
coefficients by imposing the positivity constraint with 
or without li minimization, see lH Q and IS). The 
work of [|3| extends this result to the continuous setting 
by using total variation minimization. Recently, Candes 
and Femandez-Granda showed that the solution to an 
£i -minimization problem with a — misfit 

will be close to the true signal, assuming that locations of 
any two consecutive nonzero coefficients are separated 
by at least four times the super-resolution factor n/m 
j!]. Such optimization ideas have the advantage of being 
easily generalizable to higher dimensions. On the flip 
side, li minimization super-resolution is known to fail 
on sparse signals with nearby components that alternate 
signs. 

In this paper, we discuss a simple algorithm for 
solving ([ill based on 

« subspace identification as in the matrix pencil 
method, but without the subsequent eigenvalue 
computation; and 

« a removal procedure for tightening the active set, 
remindful of a step in certain greedy pursuits. 

This algorithm can outperform the well-known matrix 
pencil method, as we show in the numerical section, and 
it is generalizable to higher dimensions. It is a one-pass 
procedure that does not suffer from slow convergence 
in situations of high coherence. We also show that 
the algorithm provides perfect recovery for the (not 
combinatoriaUy hard in the Fourier case) noiseless £o 



problem 

min|suppa;| s.t. Ax = ij. (2) 

X 

II. Noiseless subspace identification 

For completeness we start by recalling the classical 
uniqueness result for (|2). 

Lemma 1. Let Xq € R" with support T such that m > 
2\T\, and let y = Axq. Then the unique minimizer o/(|2]) 

is Xq. 

We make use of the following notations. Denote 
supp.To by T, and write At for the restriction of A to 
its columns in T. Let T"^ for the complement of T. Let 
a/c for the k-th column of A. The superscript L is used 
to denote a restriction of a matrix to its first L rows, as 
in A^. 

The "superset method" hinges on a special property 
that the partial Fourier matrix A does not share with 
arbitrary dictionaries: each column is translation- 
invariant in the sense that any restriction of at to s < m 
consecutive elements gives rise to the same sequence, 
up to an overall scalar In other words, exponentials are 
eigenfunctions of the translation operator This structure 
is important. There is an opportunity cost in ignoring 
it and treating ([T]i as a generic compressed sensing 
problem. 

A way to leverage translation invariance is to recog- 
nize that it gives access to the subspace spanned by the 
atoms a/c for k e T, such that y = J2keTi^o)k0.k- 
Algorithmically, one picks a number 1 < L < m and 
juxtaposes translated copies of (restrictions of) y into the 
Hankel matrix Y = Hankel(2/), defined as 

/ yo yi ■■■ y,n-L-i \ 

yi y2 ■■■ 2/m-L 

Y = 

\ VL-l VL ■■■ y,n ) 

The range of Y is the subspace we seek. 

Lemma 2. If L> \T\, then the rank ofY is \T\, and 

RanY = RanA!f. 

The lemma suggests a simple recovery procedure in 
the noiseless case: loop over all the candidate atoms Ofe 
for — n/2 < k < n/2 and select those for which the 
angle 

Z(a^,Rany) = 0. (3) 

Once the set T is identified, the solution is obtained by 
solving the determined system 



This procedure (unsurprisingly) provides a solution to 
the noise-free io sparse recovery problem (|2). 

Theorem 3. Let xq G K" with support T such that 
m > 2\T\, and let y ~ Axq- Consider x defined by 
di]) and (0), where the Hankel matrix Y is built with 

\T\ + l<L<m-\T\- 1. Then x = xq. 

The proofs of lemma |2] and theorem |3] hinge on the 
fact that A has full spark. 

The idea of subspace identification is at the heart of 
a different method, the matrix pencil, which seeks the 
rank-reducing numbers z of the pencil 

Y~ zY, 

where F is F with its first row removed, and Y_isY with 
its last row removed. These numbers z are computed as 
the generaUzed eigenvalues of the couple {Y_*Y,Y_*Y_). 
z can also be found via solving the eigenvalues of the 
matrix Y^Y. When |T| < L < m - |T|, the collection 
of these generalized eigenvalues includes ^'^^i^l'^ for 
A: G T, as well as m — L — \T\ zeros. There exist 
variants that consider a Toeplitz matrix instead of a 
Hankel matrix, with slightly better numerical stability 
properties. When L = \T\, the matrix pencil method 
reduces to Prony's method, a numerically inferior choice 
that should be avoided in practice if possible. 

III. Noisy subspace identification 

The problem becomes more difficult when the ob- 
servations are contaminated by noise. In this situation 
RanA^ 7^ RanF, though in low-noise situations we 
may still be able to recover T from the indices of the 
smallest angles Z(a^,Rany). 

Proposition 4. Let y = yo+e with e ~ N{0, a^Im), and 
form the corresponding L x (m — L) matrices Y and lo 
as previously. Denote the singular values ofY^ by s„,o- 
Let r = rank{Yi^)). Then there exists positive Ci,Ci and 
c, such that with probability at least 1 — cim~*^^. 



\ sin Z{aj^, Ran Yo) - sin Z{a^,RanY)\ < cei, (5) 



where 



\J rm log m 



(6) 



Atxt = y, XT" = 0. 



(4) 



In order to write bounds involving the singular values 
Si of Y instead of those of Y^, one may use Weyl's 
inequality \si — Sifl\ < \ \ Hankel(e)j|, which can in turn 
be controlled as 0(cr-y/mlogm) with high probability. 

The subspace identification step now gathers all the 
values of k compatible with the (null) hypothesis that 



Z(a^, Ran Yq) = 0, i.e., those which for some adequate 
c > obey 

sinZ(a^,RanF) < C£i. 

The resuhing set Vl of indices is only expected to be a 
superset of the true support T, with high probabiUty. 

A second step is now needed to prune Vl in order 
to extract T. For this purpose, a loop over k is set up 
where we test the membership of y in RanAn^j,, the 
range of Aq_ with the fc-th column removed. We are now 
considering a new set of angles where the roles of y and 
A are reversed: in a noiseless situation, fc S T if and 
only if 

Z(2;,RanAf^\fc)^0. 

When noise is present, we first filter out the noise off 
by projecting y onto the range of A^, then estimate 
k E T only when the angle is above a certain threshold. 
It is easier to work directly with projections 11: 

llnny - nn\fc?;|| = sinZ(no?/,RanAo\fe) \\lliw\\- 

The effect of noise on the left-hand side is as follows. 

Proposition 5. Let y = yo + e with e ^ 7V(0,cT^/m). 
Let Tiny be the projection of y onto Ran An, and let 
An = Tin 

with high probability. 



Algorithm 1 Superset selection and pruning 



na\fc. Then there exists c > such that. 



|||Any||-|lAnyo|l|<c£2, 

with £2 = cr. 

Algorithm [1] for the superset method implements the 
removal step in an iterative fashion, one atom at a time. 

IV. Experimental Results 

In the first simulation, we fix n = 1000 and m = 120 
and construct an n-dimensional signal Xq whose nonzero 
components are well separated by at least An/m, a 
distance equivalent to four times the super-resolution 
factor n/m. The spike magnitudes are independently set 
to ±1/V29 with probability 1/2. The noise vector e 
is drawn from N{0,a^lm) with a = lO^'^. We fix the 
thresholds ei via ^ with c = 1 and £2 = lOcr. Through- 
out our simulations, we set L = [?7i/3j. As can be 
seen from Fig. [T] top row, the recovered signal from the 
superset method is reasonable, with ||a; — a:^o|l2 — 0.075, 
while the reconstruction via ^1 -minimization tends to 
exhibit incorrect clusters around the true spikes. 

Our next simulation considers a more challenging 
signal model with a strongly coherent matrix A. For 
example, with n = 1000 and = 120, the coher- 
ence of the matrix A with normalized columns a; is 
fi = max,;^j I (a,;, flj) I = 0.9765. The signal in this 



input: Partial Fourier matrix A G (7™^", 
e, parameter L, thresholds £1 and £2- 
initialization: Y = Hankel(y) e C^^f™" 
support identification 
decompose: QR = YE, Q e 
project: a^- A^j^y ( for all k) 



y = Axq + 



••Lxr 



Ik 



ak 



'ak 



/hk\ 



while true do 
decompose: 

remove: 



= {fc : 7fc < £1} 
QR = AnE, Q e C™^l^^l 

QQ*)y\\2 



Vfc e n-. Q^k)R(k) 
Sk ^ \\{Q(k)Qlk) 

fco ^ argminj, Sk 

if 4o <e2,fl^ fl\k„ 

else break 



^n\kE{k) 



end while 
output: X 



argmm^ \\y ■ 



AnxW 



Our proposed algorithm 



LI minimization 



— ' True signal 

— Recovered signal 



□0 400 600 f 
Signai dimension 
Our proposed algorithm 



— ^ True signal 
Recovered signai 



Signai dimension 
LI minimization 



It 


— ^ True signai 

— Recovered signal 


I 











— ^ True signal 

— Recovered signal 




— p 



Fig. L Original (blue) and recovered (red) signals. Left 
column: the superset method. Right column: ^i-minimization. 
Top row: a signal with well-separated spikes. Bottom row: 
spike spacing below the Rayleigh length. 



simulation is shown in Fig. [T] bottom row. It consists of 
five spike clusters: each of the first two clusters consists 
of a single spike, and each of the last four clusters 
contains two neighboring spikes. The signs of these 
neighboring spikes either agree or differ We set to, a and 
£2 as in the previous simulation, and we let the constant 
c in the equation ^ of £1 equal to 5. Recovery via the 
superset method is accurate, while £1 minimization fails 
at least with clusters of opposite-sign spikes. 

In the next simulation, we consider a signal of 
size n ~ 1000 which contains two nearby spikes 
at locations [100,101] and has magnitudes 1/V^ and 



— 1/\/2. We empirically investigate the algorithm's abil- 
ity to recover the signal from varying measurements 
m = {10, 20, 220} and noise levels logioa = 
{— 3.5, — 3.4, — 2}. For each pair (m,(T), we report 
the frequency of success over 100 random realizations 
of e. The greyscale goes from white (100 successes) to 
black (100 failures). A trial is declared successful if the 
recovered x satisfies ||iz; — a-pHj / ||a;o|l2 < 10"'^. The 
horizontal axis indicates the noise level a in log scale, 
and the vertical axis indicates logj^Q(l — /i) where i^l is 
the coherence as earlier 

We note that the coherence is inversely proportional to 
the amount of measurements m and proportional to the 
super-resolution factor n/m: increasing m (decreasing 
the super-resolution factor) will reduce the coherence fx. 
On the vertical axis, smaller values imply higher coher- 
ence, or equivalently smaller amount of measurements. 
As shown in Fig. |2] for reasonably small noise, the 
algorithm is able to recover the signal exactly even the 
coherence is nearly 1. 

For reference, we also compare the superset method 
with the matrix pencil method as set up in ifTOl . The 
noise is filtered out by preparing low-rank approxima- 
tions of Y_ and Y where only the singular values above 
ca\JL logi are kept, for some heuristically optimized 
constant c. Two more signals are considered: (1) a 3- 
spase signal consisting of three neighboring spikes, each 
of magnitude 1 / \/3 with alternating signs, and (2) a 4- 
sparse signal with neighboring spikes of alternating signs 
and equal magnitude 1/2. Fig.|2]is a good illustration of 
the contrasting numerical behaviors of the two methods: 
the matrix pencil is often the better method in the special 
case of a signal with 2 spikes, but loses ground to the 
superset method in various cases of progressively less 
sparse signals. Understanding the performance of the 
matrix pencil would require formulating a lower bound 
on the (typically extremely small) 5-th eigenvalues of 
Yb where S is the sparsity of yo- 

V. Conclusion 

Empirical evidence is presented for the potential of 
the superset method as a viable computational method 
for super-resolution. Further theoretical justifications will 
be presented elsewhere. 
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